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Abstract 

The Osher-Chakrabarthy family of linear flux-modification schemes is considered. 
Improved lower bounds on the compression factors are provided, which suggest the 
viability of using the unlimited version. The LLF flux formula is combined with 
these schemes in order to obtain efficient finite-difference algorithms. The resulting 
schemes are applied to a battery of numerical tests, going from advection and Burg- 
ers equations to Euler and MHD equations, including the double Mach reflection and 
the Orszag-Tang 2D vortex problem. Total-variation-bounded behavior is evident in 
all cases, even with time-independent upper bounds. The proposed schemes, how- 
ever, do not deal properly with compound shocks, arising from non-convex fluxes, 
as shown by Buckley-Leverett test simulations. 



Key words: Hyperbolic conservation laws. Numerical methods 
PACS: 47.11Df, 47.11BC 



Work jointly supported by the European Commission FEDER funds, the Spanish 
Ministry of Science and the Balearic Islands Government 

^ C. Bona acknowledges the Charles University in Prague for his hospitality during 
the completion of this work and specially Tomas Ledvinka for useful suggestions 
and discussions. 

^ C. Bona-Casas acknowledges the support from the FPU/2006-02226 fellowship of 
the Spanish Ministry of Science and Education 

^ J. Terradas acknowledges the support from the Research Council fellowship 
F/06/65 of the Katholieke Universiteit Leuven 



Preprint submitted to Elsevier 



12 December 2008 



1 Introduction 



The study of hyperbolic conservation laws, as described by 

dtu + dj{u) = 0, (1) 

is a classical topic in Computational Fluid Dynamics (CFD). We have noted 
here by m a generic array of dynamical fields, and we will assume strong 
hyperbolicity, so that the characteristic matrix 

A{u) = df/du (2) 



has real eigenvalues and a full set of eigenvectors. 

As it is well known, the system ([1]) admits weak solutions, so that the com- 
ponents of u may show piecewise-smooth profiles. Standard finite-difference 
schemes, like the Lax-Wendroff [Tj or MacCormack [2] ones, produce spurious 
overshots and oscillations at non-smooth points which can mask the physical 
solutions, even leading to code crashing. These deviations do not diminish 
with resolution, in analogy with the Gibbs phenomenon found in the Fourier 
series development of discontinuous functions. 

This difficulty was overcome in the pioneering work of Godunov [3]. On a 
uniform computational grid Xj = jAx, equation ([T]) can be approximated by 
the semi-discrete equation 

dtuj = --^ {hj+i/2 - /ij-1/2) , (3) 



where the interface flux /ij+1/2 is computed by an upwind-biased formula from 
the neighbor grid nodes. In the scalar case, one can define the total variation 
of a discrete function as 

TViu) = Y: • (4) 

j 



In the case of systems, the total variation is defined as the sum of the total 
variation of the components. Godunov scheme is total-variation-diminishing 
(TVD), meaning that TV{u) does not increase during numerical evolution. It 
is obvious that TVD schemes can not develop spurious oscillations: monotonic 
initial data preserve their monotonicity during time evolution. Moreover, the 
TVD property can be seen as a strong form of stability: any blow-up of the 
numerical solution is excluded, as far as it would increase the total variation. 
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Godunov scheme is the prototype of the so-called upwind-biased schemes, 
which require either the exact or some approximate form of spectral decom- 
position of the characteristic matrix ([2]). This makes them both computation- 
ally expensive and difficult to extend to the multidimensional case. A much 
simpler alternative is provided by the local Lax-Friedrichs (LLF) scheme (Ru- 
sanov scheme 

hj+l/2 = ^ [fj+l + fj - \j+l/2 {Uj+1 - Uj) ] , (5) 

where A is the spectral radius of the characteristic matrix and we have taken 
Aj+i/2 = max{Xj, Xj+i) . (6) 

The LLF scheme is the prototype of the so-called centered schemes, which 
deserve a revived interest nowadays in view of multidimensional applications. 
It clear from ([3l [5]) that the LLF scheme, like the Godunov one, is only first- 
order accurate in space. Second-order accuracy can be obtained following the 
Harten modified-fiux approach [5], which was soon extended to very-high ac- 
curacy (up to 15th order) by Osher and Chakrabarthy [6]. The basic idea is 
to replace the lower order TVD fiux /ij+1/2 by a modified fiux fj+1/2, obtained 
by some interpolation procedure involving a higher number of nodes. 

All these high-resolution schemes require some form of fiux-correction limiters 
in order to ensure the TVD property. As a consequence, accuracy is reduced to 
(at most) first order at non-sonic critical points, where the limiters come into 
play. In order to circumvent this problem, one can relax the TVD condition, 
demanding instead that the total variation is bounded, that is 

TV{u) <B , (7) 

where the upper bound B is independent of the resolution, but could depend 
on the elapsed time. Even if we are ready to relax the stronger TVD require- 
ment, keeping the bound ([7]) is important from the theoretical point of view. 
One major advantage of total-variation-bounded (TVB) schemes is that there 
is a convergent (in Lj^^) subsequence as Ax — > to a weak solution of ([1]). 
If an additional entropy condition is satisfied, then the proposed scheme is 
convergent (see for instance Ref. ^). 

An interesting example of such TVB schemes was given by Shu [S] , by softening 
the flux limiters proposed in Ref. [6] . Although the TVB property is proven for 
the schemes presented in [8], based on a linear flux- modification procedure, 
a rigorous proof is still unavailable for more complex cases. An important 
example is provided by the essentially-non-oscillatory (ENO) methods [H] PH] . 
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where the TVD property is relaxed in a different way. Numerical evidence 
shows that ENO schemes, as well as their weighted-ENO variants p!T]-|13j. 
deserve their name: the TVB property is satisfied in practice, even with time- 
independent bounds. An implementation of these high-resolution methods for 
the LLP Flux is given in Refs. IH] |15]. 

In this paper we consider (the unlimited version of) the '/5-family' of Osher- 
Chakrabarthy linear flux-modification algorithms in the semidiscrete frame- 
work. The choices of the (3 parameter are optimized for their use with the 
unlimited version, by refining the bounds given in the original paper [6]. As 
a lower order TVD flux, we will use the standard LLP one ([5]). The resulting 
high-resolution schemes can then be recast into a compact finite-difference for- 
mula. This greatly increases computational efficiency, specially with a view to 
three-dimensional applications. We perform a battery of standard tests in one 
space dimension, covering advection. Burgers and Euler equations, in order 
to show that the TVB property is fulfilled in practice for the selected values 
of the (3 parameter. We are not able, however, of getting the right result for 
compound shocks, arising from non-convex fluxes; this is illustrated by the 
Buckley-Leverett test simulations. We also consider some multidimensional 
tests cases with the Euler and magneto-hydrodynamics (MHD) equations, in- 
cluding the double Mach reflection and the Orszag-Tang 2D vortex problem. 
Total-variation-bounded behavior is evident in all the proposed cases, even 
with time-independent upper bounds. 



2 The Osher-Chakravarthy /3-schemes 

PoUowing Ref. fB], let us consider the centered 2m — 1 order schemes: 

d,u, = -C^^f, + {-ir-^(3{^xf^-'D^D^-\dfl_,/, - dfr^,,,) , (8) 



where C is the central 2mth-order-accurate difference operator with a stencil 
of 2m + 1 grid points, and we have used the standard notation D± for the 
elementary difference operators. The flux differences df^ are defined as follows: 

= fj+^ ~ ^/j>i/2 = ^i+i/2 - fj > (9) 

where is any (lowest-order) TVD flux. The (3 parameter in the dissipa- 

tive term in ^ is assumed to be positive: a necessary condition for stabihty. 
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The algorithms ([8]) can be put into an exphcit flux-conservative form by re- 
placing the lowest order flux /ij+1/2 in © by [6] [8] 



m— 1 



/i+1/2 - ^i+1/2 + {'^"k dfj+k+1/2 + ^"l df. 



k=—m+l 



j+k+1/2 



(10) 



where 



d1 = y'l - {-^ffi 
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The TVD property is enforced by limiting the flux differences df^ (see [6], [8] 
for the details). As a generic example, (i/jt|_fc+i/2 in din]) is replaced by 



minmod{ c^/^^fc+i/a, M/^>i/2> 



(15) 



where 6 is a compression factor. This replacement introduces a non-linear 
component in the linear flux-correction formula (llUp . The resulting scheme 
will be TVD if and only if: 



"sr~^ m. dfj+k+1/2 df-^j^_^i2 



k= — m+l 
m—1 



df,+. 



> 



i+1/2 



n - 1 M \^ ^ '^/i+fc+1/2 ^/j+fc-1/2 ^ 

k=-m+l •'i-1/2 



(16) 

(17) 
(18) 



where we have assumed a time discretization based on the forward Euler step, 
so that the last condition provides an upper bound on the time step At. 
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3 Compression factor optimization 

In the original paper [6], the ansatz 



2m , 
' m 



(19) 



was used for getting a sufficient condition from (fT6l[T7I) . amounting to a simple 
constraint on the range of the compression parameter b 



< 6 < [ 1 + 2/3 




m 1 



1-1 



(20) 



Allowing for ( l20i) . the upper bound hmax increases with /?, which is in turn 
bounded by (fT9|) . For the third-order scheme (m = 2), the optimal choice 
would then be /3 = 1/12, so that the compression parameter may reach hmax = 
4, still preserving the TVD property. This means that, for monotonic profiles, 
the fiux-correction limiters would act only where the higher order corrections 
in neighboring computational cells differ at least by a factor of four. This is not 
to be expected in practical, good resolution, simulations of smooth profiles, 
even when large gradients appear. This high-compression-factor property can 
be at the origin of the robust behavior of these schemes, even in their unlimited 
form, as we will see in the numerical applications presented below. 

As far as we are proposing to use the unlimited version, it makes sense to 
find the choices of /? that maximize the compression factor, going beyond 
the ansatz ( fT9|) . Higher values of hmax can be actually obtained by a detailed 
case-by-case study of the original TVD conditions (fT6| [T7|) . For instance, by 
reordering the terms in (fT7|) we get 

Dj-i/2 = 1 + E (4-1 - dn -jS^ > , (21) 



where we assume djj^ = when \k\ > m. A sufficient condition for ( 1211) to 
hold is 

1 + - rfo" + & E min{dn - 4", 0) > , (22) 



6 



which actually refines the former condition (120|) . The same reasoning shows 
that, allowing for (fTTj) . a sufficient condition for (ITSl) to hold is: 

A, ^ [ rf-'T - rfo'" + & E max{d,"^, - d,^, 0) ] < 1/2 . (23) 



For the simpler non-trivial cases we have (decreasing k order): 



d,' = {/3--^, i-2/3, /3 + 1) (24) 



For m = 2, condition (122|) leads then to: 



b < 7/2 + 18/3 W<J^) (26) 

, < (i<^<^). (27) 

It follows that the optimal values for the third-order scheme are 

^ = ^ , bmax = 5 . (28) 



For the fifth-order scheme (m = 3), condition fl22l) leads instead to: 



37 + 600/3 1^ 

< 37 + 600^ (1</,<1) (30) 

- 15 + 60/3 ^60-^-75^ ^ ^ 

37 + 600/3 ,2 ^ 37, 

6 < (— < 8 < ). (31) 

- 7 + 360/3 ^5-^- 600^ ^ ' 

It follows that the optimal values for the fifth-order scheme are 



Note that the ansatz (1191) gives a smaller compression factor bmax = 9/4 and, 
more important, the optimal (3 value in this case is beyond the original bound 
1/60. Note also that the values of the compression parameter tend to diminish 
with the accuracy order of the algorithm. This suggests that higher-order cases 
m > 3 may not be so useful in the unlimited case. 
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4 Finite difference version 



The linear flux-modification scheme described in the preceding section can be 
applied to any lower-order TVD flux. The case of the LLF flux has actually 
been considered in [S]. Our objective here is to obtain a scheme which can be 
cast as a simple finite-difference algorithm, so that we will take advantage of 
the simplicity of the LLF flux ([5]), which can be written in flux- vector-splitting 
(FVS) form: 

hj+i/2 = // + f~+i , ft = \ ifi ± ^i±i/2 Uj] . (33) 

The FVS form ( l33l) . like the original one ([5]), is just first-order accurate. We 
will extend it to higher-order accuracy by means of the Osher-Chakrabarthy 
algorithm, as described in the previous sections. The flux differences (|9]) in 
this case get the simple form: 

= 1/2 [ fj+i - fj ± Aj+i/2 K+i -Uj)] . (34) 

The linear character of this formula allows to get a compact finite-difference 
expression for the whole scheme. Allowing for the semi-discrete algorithm 
([8]) can be written as 

dtu, = -C^^^f^ + (-l)"^-i/3(Ax)2™-iZ};"Z):['-i(A,_i/2 D_u,) , (35) 

which amounts to assume a 2mth-order-accurate central difference operator 
acting on the flux terms plus a dissipation operator of order 2m depending 
on the spectral radius A. As we will see below, the resulting finite-difference 
scheme ( 135|) provides a cost-effective alternative for CFD simulations. 

Let us remark here that the choices (l28l [32]) derived in the previous section are 
optimal for a generic choice of the lowest-order TVD Flux. In the LLF case 
([S]), however, it is clear that the spectral radius can be multiplied by a global 
magnifying factor K > 1, while keeping the TVD properties. Allowing for the 
finite-difference form fl35|) of the unlimited version, magnifying A amounts to 
magnify /5, that is: 

(/3, KX) ^ {K/3, A) . (36) 

It follows that the values of the compression factor bmax obtained in the previ- 
ous section must be interpreted just as lower-bound estimates. In particular, 
the equivalence (l36l) implies that any compression factor bound obtained for 
a particular value j3o applies as well to all values (3 > (3q. This agrees with 
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the interpretation of the second term in (l35l) as modelhng numerical dissipa- 
tion. On the other side, this dissipation term is actually introducing the main 
truncation error. We will use then in what follows the (3 values in (1251 [5^ . 
which are still optimal in the sense that they provide the lower numerical error 
compatible with the highest lower-bound for the compression parameter. 



5 One-dimensional numerical tests 

We will test now the behavior of the centered finite-difference scheme f l35l) 
by means of some standard numerical experiments in one space dimension. 
We will use here the well-known method-of-lines (MoL) fT6] in order to deal 
separately with the space and the time discretization. In every case, the time 
discretization will be implemented by the following strong-stability-preserving 
(SSP), 3rd-order-accurate, Runge-Kutta algorithm [T7] : 

u* = E{u'') 

u** = ^u'' + - Eiu*) (37) 

where E{u) is the basic Euler step, that is 

E{uj) = Uj + /\t {dtUj) , (38) 

and dt Uj is computed by the finite difference formula fl35l) . 
5.1 Advection equation 

Let us start by the scalar advection equation. This is the simplest linear case, 
but it allows to test the propagation of arbitrary initial profiles, containing 
jump discontinuities and corner points, departing from smoothness in many 
different ways. This is the case of the Balsara-Shu profile [11], which will be 
evolved with periodic boundary conditions. 

We compare in Fig.[T]the numerical result with the exact solution after a single 
round trip, for two different resolutions. The third-order five-points formula 
from the proposed class (135!) has been used with /3 = 1/12 in both cases. The 
propagation speed in the simulation agrees with the exact one, as expected 
for a third-order-accurate algorithm. The smooth regions are described cor- 
rectly: even the height of the two regular maxima is not reduced too much 
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Fig. 1. Advection of the Balsara-Shu profile in a numerical mesh of either 400 points 
(upper panel) or 800 points (lower panel). A third-order scheme (m = 2, /? = 1/12) 
is used in both cases. The results are compared with the initial profile (dotted line) 
after a single round-trip. 

by dissipation, as expected for an unlimited algorithm with just fourth-order 
dissipation. There is a slight smearing of the jump slopes, as usual for contact 
discontinuities, which gets smaller with higher resolution. 




Fig. 2. Same as in Fig. [H but using a fifth-order scheme (m = 3) with [5 = 2/75 
(upper panel). In the lower panel we show the results after ten round-trips. The 
same settings are used in both cases. 
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Fig. 3. Advection equation. Time evolution of the total variation. The horizontal 
axis corresponds to the exact solution: TV{u) = 8. From top to bottom: m = 3 
scheme with 400 points, m = 2 scheme with 800 points, and m = 2 scheme with 
400 points. After the initial increase, which depends on the selected method, the 
TV tends to diminish. Increasing resolution just reduces the TV diminishing rate. 

Concerning monotonicity, it is clear that the total variation of the initial profile 
has increased by the riddles besides the corner points and, more visibly, near 
the jump discontinuities. By comparing the two resolutions, we see that the 
height of the overshots does not change. This means that, as in the case of the 
Gibbs phenomenon, there is no convergence by the maximum norm, although 
convergence by the L2 or similar norms is apparent from the results. On the 
other hand, it is clear that the total variation is bounded for this fixed time, 
independently of the space resolution or, equivalently, the time step size. This 
is precisely the requirement for TVB. 

We show in Fig. [2] the same simulation, in a 400 points mesh, for the fifth- 
order method (m = 3, /3 = 2/75). In the upper panel, corresponding to a 
single round-trip, we can see that one additional riddle appears at every side 
of the critical points, due to the larger (seven point) stencil. We show also in 
the lower panel the results of the same simulation after ten round-trips. The 
cumulative effect of numerical dissipation is clearly visible: the extra riddles 
tend to diminish. The total variation is not higher than the one after a single 
round trip. This statement can be verified by plotting, as we do in Fig. [HI the 
time evolution of TV{u) for the different cases considered here. In all cases, 
a sudden initial increase is followed by a clear diminishing pattern. These 
numerical results indicate that the bound on the total variation is actually 
time-independent, beyond the weaker TVB requirement. 



5.2 Burgers equation 

Burgers equation provides a simple example of a genuinely non-linear scalar 
equation. A true shock develops from smooth initial data. We will compute 
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Fig. 4. Burgers equation: evolution of an initial sinus profile. The numerical solu- 
tion (point values) is plotted versus the exact solution (continuous line), for 100 
points and 200 points resolution (left and right panels, respectively) and for the 
(m = 2, /? = 1/12) and the (m = 3, /? = 2/75) schemes (upper and lower panels, 
respectively) . 

here the evolution of an initial sinus profile, with fixed boundary conditions. 
We plot in Fig. Hlthe numerical solution values versus (the principal branch of) 
the exact solution, at the time where the shock has fully developed. We com- 
pare 100 points with 200 points resolution (left and right panels, respectively), 
and also the 3rd-order and 5th-order schemes described previously (upper and 
lower panels, respectively). Concerning the resolution effect, we can see here 
again that the spurious oscillations affect mainly the points directly connected 
with the shock, in a number depending on the stencil size but independent of 
the resolution. 

These conclusions are fully confirmed by a second simulation, obtained by 
adding a constant term to the previous initial profile, that is 



with periodic boundary conditions. We can see in Fig. [5] that a shock again 
develops, but it does no longer stand fixed: it propagates to the right. Note 
that the plot shown corresponds to t = 7. We can confirm in this case that 
both the number of spurious riddles and the magnitude of the overshots do 




(39) 
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Fig. 5. Same as in the previous figure, but now for a moving sinus profile. The 
numerical solution (point values) is plotted versus the exact solution (continuous 
line), for 100 points and 200 points resolution (left and right panels, respectively) 
and for the (m = 2, /? = 1/12) and the (m = 3, /5 = 2/75) schemes (upper and 
lower panels, respectively). 

not increase with resolution, although it is larger in this case than in the static 
shock one. We can confirm also that these effects increase with the order-of- 
accuracy of the scheme: the larger stencil adds one more riddle at every side 
and slightly larger overshots. 



These results clearly indicate convergence in the Li or similar norms (but of 
course not in the maximum norm) . Let us actually perform a convergence test 
by considering the initial profile [18] 



u{x^ 0) = 1 + - sm(7r x) 



(40) 



which is smooth up to t = 2/7r. We show in Tabled] the errors at time t = 0.3, 
where the shock has not yet appeared. The first group of values corresponds 
to the third-order method, and this is confirmed by the data both in the Li 
and the L^o norms. The second group of values corresponds to the fifth-order 
method, but only third-order accuracy is obtained from the numerical values. 
This is because we keep using the third-order Runge-Kutta algorithm (1371) for 
the time evolution. In order to properly check the space discretization accu- 
racy, we include a third group of values, obtained with the same algorithm, but 
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with a much smaller time step in order to lower the time discretization error: 
the leading error term is then due to the space discretization and the expected 
fifth order accuracy is confirmed by the numerical results, although the 
norm shows a slightly decreasing convergence rate for the higher resolution 
results. 



Nx 


Li error 


Li order 


Loo error 


Loo order 


160 


7.22579 E-6 


2.998 


5.17334 E-5 


2.981 


320 


9.04719 E-7 


2.999 


6.55306 E-6 


2.994 


640 


1.13182 E-7 


3.000 


8.22735 E-7 


2.998 


1280 


1.41486 E-8 




1.03006 E-7 




160 


1.44981 E-6 


3.017 


9.57814 E-6 


2.981 


320 


1.79043 E-7 


3.005 


1.21318 E-6 


2.997 


640 


2.23035 E-8 


3.003 


1.51957 E-7 


2.999 


1280 


2.78216 E-9 




1.90041 E-8 




160 


7.09726 E-8 


4.88 


8.6567 E-7 


4.97 


320 


2.41410 E-9 


4.76 


2.76804 E-8 


3.98 


640 


8.92936 E-11 


4.91 


1.75192 E-9 


3.48 


1280 


2.95859 E-12 




1.36890 E-11 





Table 1 

Burgers problem. Norm of the errors and convergence rate at t = 0.3 for the initial 
profile ()40p . The first group of values corresponds to the m = 2 method with At = 
O.6A3;. The second group corresponds to the m = 3 method with the same time step. 
The third group corresponds again to the m = 3 method, but with At = 0.06Ax. 



5.3 Buckley - Lev erett problem 



A more demanding test, still for the scalar case, is provided by the Buckley- 
Leverett equation which models two-phase fiows that arise in oil-recovery prob- 
lems [7]. This equation contains a non-convex (s-shaped) fiux of the form 



Non-convex fiuxes can lead to compound shock waves which are shocks adja- 
cent to a rarefaction wave with wave speed equal to the shock speed at the 
point of attachment. 
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Fig. 6. Buckley-Leverett's problem. The continuous line corresponds to the LLF 
first-order algorithm, with 10.000 points, as a replacement for the exact solution. 
The crosses line corresponds to the third-order algorithm (m = 2, f3 = 1/12) with 
200 points, converging towards a different solution. 

We will perform first a simulation with the initial data 

,0 < X < 1 - l/y/2 

u{x) = { - (42) 

1 1 - 1/V2 < X <1 



so that the inflexion point in the flux (HTj) lies inside the interval spanned by 
the data. 

The exact solution in this case is well approximated by a very-high-resolution 
(10.000 points) simulation using the flrst-order LLF algorithm, as displayed 
in Fig. [6] (continuous line). We see a right-propagating compound shock wave, 
consisting of a shock followed by a rarefaction wave, which propagates in the 
same direction. The results for our third-order algorithm, represented by the 
crosses line in Fig. El fail to reproduce correctly the rarefaction wave, which 




Fig. 7. Same as in the previous figure, but now for two different dynamical ranges, 
which avoid the flux inflexion point. In the left panel, an ordinary rarefaction wave 
appears, which is correctly modelled by the third-order algorithm. In the right panel, 
a simple shock appears, well captured by the third-order algorithm. 
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is replaced by an spurious intermediate state, resulting into a slower shock 
propagation speed. 

In order to single out the problem, we have performed simulations for the 
same flux fHT]) but with a dynamical range that avoids the inflexion point 
either from below or from above. The results are plotted in Fig. [TJ where we see 
either an ordinary rarefaction wave (left panel) or a simple shock (right panel), 
but no compound shock. In both cases, the third-order algorithm (m = 2, 
/3 = 12) is able to model correctly the dynamics. This results indicate that the 
problem with compound shocks can be triggered by the presence of overshots 
at the connection point between the shock and the associated rarefaction wave, 
which can break the compound structure. The TVD character of the LLF flux 
prevents this problem to arise, as it is clearly shown in Fig. |6] (continuous line). 



5.4 Euler equations 



Euler equations for fluid dynamics are a convenient arena for testing the pro- 
posed schemes beyond the scalar case. In the ideal gas case, we can check the 
numerical results against well-known exact solutions containing shocks, con- 
tact discontinuities and rarefaction waves. We will deal first with the classical 





Fig. 8. Sod shock tube problem. Density and speed profiles (left and right panels, 
respectively), for the (m = 2, /3 = 1/12) and the (m = 3, /3 = 2/75) schemes (upper 
and lower panels, respectively). 
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Sod shock-tube test [19] with a standard 200 points resolution. 



We plot in Fig. [8] the gas density and speed profiles (left and right panels, 
respectively). Looking at the 3rd-order scheme results (upper panels), we see 
that both the rarefaction wave and the shock are perfectly resolved, whereas 
the contact discontinuity is smeared out. As a consequence, the main overshots 
are just besides the shock, specially visible in the speed profile, where the jump 
is much higher. Concerning the 5th-order scheme (lower panels), the contact 
discontinuity is slightly better resolved. This is however at the price of extra 
riddles and more visible overshots, so that the 3rd-order scheme seems to be 
more convenient. 

A more demanding test is obtained when assuming a discontinuity in the initial 
speed, as in the Lax test [20j. As we see in Fig. [9l we get the same behavior 
than for the Sod test case. The main difference is that the density jump at 
the contact discontinuity is much higher: the smearing of the density profile 
there is more visible, in contrast with the sharp shock profile nearby. Note 
also that some speed overshots are greater than the ones arising in the Sod 
test case (we have kept here the same 200 points resolution for comparison). 
The third-order algorithm seems to be more convenient again in this case. 
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Fig. 9. Lax shock tube problem. Density and speed profiles (left and right panels, 
respectively), for the (m = 2, /? = 1/12) and the (m = 3, /? = 2/75) schemes (upper 
and lower panels, respectively). 
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6 Multidimensional tests 



The results of this paper can be extended to a multidimensional case in a 
simple way. The semi-discrete equation ([3]) can be written in a rectangular 
grid as follows: 

dtUiJ = (/i+l/2,j - fi-l/2,j) - ifiJ+l/2 - fi,j-l/2), (43) 



and the numerical flux can be computed by applying (fTOl) to every single 
direction. Note however that the restriction ( l23l) on the time step must be 
extended in this case to 

A, At + d - - d - + 6 ^ max(dA - 4", 0) ] < 1/2 . (44) 



In the finite-difference version (135|1 . the extension to the multidimensional 
case amounts to replicate the right-hand-side difference operators for every 
single direction: no cross-derivative terms are required. This multidimensional 
extension allows to deal with some MHD tests, which add more complexity 
to the dynamics, clearly beyond the simple tests considered in the previous 
section. 



6.1 The Orszag-Tang 2D vortex problem 



As a first multi-dimensional example, let us consider here the Orszag-Tang vor- 
tex problem pTj . This is a well-known model problem for testing the transition 
to supersonic magnetohydrodynamical (MHD) turbulence and has become a 
common test of numerical MHD codes in two dimensions. 

A barotropic fluid (7 = 5/3) is considered in a doubly periodic domain 
[0, 27r]^, with uniform density p and pressure p. A velocity vortex given by 
V = (— siny, sinx), corresponding to a Mach 1 rotation cell, is superimposed 
with a magnetic field B = (— siny, sin2x), describing magnetic islands with 
half the horizontal wavelength of the velocity roll. As a result, the magnetic 
field and the flow velocity differ in their modal structures along one spatial 
direction. 

In Fig. [To] (left panel) the temperature, T = pj p, is represented at a given time 
instant (t = 3.14). The figure clearly shows how the dynamics is an intricate 
interplay of shock formation and collision. The numerical scheme, with m = 2 
and /? = 1/12 seems to handle the Orszag-Tang problem quite well. In Fig. [Tn] 
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Fig. 10. Temperature at t = 3.14 in the Orszag-Tang vortex test problem. In this 
simulation the grid has 200 x 200 mesh points. In the left panel the third-order 
scheme (m = 2, f3 = 1/12) has been used while in the right panel the result is for a 
a second order scheme built from the Roe- type solver and the MC limiter. 

(right panel) we plot the results for the same problem using a second order 
scheme built from the Roe-type solver and the monotonized-central (MC) 
symmetric limiter [22]. The results with both methods are qualitatively very 
similar. 



6.2 Torrilhon MHD shock tube problem 



We now consider the MHD shock tube problem described by Torrilhon [23] 
to investigate dynamical situations close to critical solutions. We will assume 
again a barotropic fluid with 7 = 5/3. The initial conditions for the compo- 
nents of the magnetic field (52,-83) are (cos ^, sin ^), with 6* = for x < 0. 
Depending on the angle 6 between the left and right transverse components of 
the magnetic field, different types of solutions are found. Regular r-solutions 
consist only of shocks or contact discontinuities. Critical c-solutions appear 
in the coplanar case, where the angle 9 is an integer multiple of vr. These 
solutions can contain also non-regular waves, such as compound waves. 

We consider the situation for an almost co-planar case, ^ = 3. Analytically, 
this has a regular r-solution, but the numerical solution is attracted towards 
the nearby critical solution for 6 = n. Fig. [TT] shows the density profile plotted 
together with the correct r-solution (solid black line) and the co-planar c- 
solution (dashed line). The r-solution has, from left to right, a rarefaction, a 
rotation, a shock, a contact discontinuity, a shock, a rotation and a rarefaction. 
The discrepancies among the different solutions are mainly in the interval 
[-0.35,-0.1]. 
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Fig. 11. Plot of the density p aXt = 0.4 for the almost co-planar problem with = 3. 
In this simulation 5000 mesh points have been used. The dashed line represents the 
critical c-solution while the solid black line is the correct r-solution. Both solutions 
differ clearly in the interval [—0.35,-0.1]. The numerical simulation lies between 
the two. 



0.90 




Fig. 12. Same as Fig. llll but enlarging the interval where the discrepancies show up. 
In addition to the exact regular and critical solutions, we plot, from top to bottom, 
the simulations for schemes using LLF with minmod limiter, LLF with MC limiter, 
the unlimited m = 2 algorithm, a Roe solver with MC limiter and the unlimited 
m = 3 algorithm. 

This interval is magnified in Fig. [T21 The solid black line is the correct r- 
solution while the dashed line represents the critical c-solution. We see that 
the solutions with m = 2 and m = 3 tend to the correct solution although 
they keep some remnant from the c-solution. For comparison purposes we have 
also represented the numerical solution obtained with other schemes. We have 
used a second order LLF scheme and a second order Roe solver with either 



20 



the minmod or the MC slope hmiters. 

The LLF scheme with the minmod hmiter gets too close to the c-solution, 
even for this high-resolution simulation. The situation improves by replacing 
the minmod limiter by the MC one, but still gets farther from the right solution 
than the schemes proposed in this paper. Only the combination of a Roe-type 
solver with the MC limiter improves the results of the third-order scheme 
(m = 2), but not those of the fifth-order scheme (m = 3). This problem 
provides one specific example in which the fifth-order scheme seems to be more 
convenient than the third order one: the extra riddles are actually compensated 
by a clear improvement in the solution accuracy. 

6.3 Double Mach reflection problem 

This problem is a standard test case for high-resolution schemes. It corre- 
sponds to an experimental setting in which a shock is driven down a tube which 
contains a wedge. We will adopt here the standard configuration proposed by 
Woodward and Colella [21], involving a Mach 10 shock in air (7 = 1.4) at 
a 60° angle with a reflecting wall. The air ahead of the shock is stationary 
with a density of 1.4 and a pressure of 1. The reflecting wall lies at the bot- 
tom of the computational domain, starting at a; = 1/6. Allowing for this, the 
exact post-shock condition is imposed at the bottom boundary in the region 
< a; < 1/6 and a reflecting wall condition is imposed for the rest. Inflow 
(post-shock) conditions are used at the left and top boundaries, whereas an 
outflow (no gradient) condition is used for the right boundary. 




o 0.5 1 1.5 2 2.5 3 



Fig. 13. Double Mach reflection. Density plot at t = 0.2. The simulation is made 
with the 3rd-order method (m = 2, /3 = 1/12) with Ax = Ay = 1/240. 30 evenly 
spaced density contours are shown. 
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Fig. 14. Same as Fig. [T3l but enlarging the lower right corner. The top panels 
correspond to the third-order method (m = 2, /? = 1/12), with resolution of either 
1/240 (left) or 1/480 (right). The bottom panels show the same for the 5th-order 
method (m = 3, (3 = 2/75). Both the jet near the bottom wall and the weak 
shock, generated at the kink in the main reflected shock, are well resolved. A vortex 
structure at the bottom of the diagonal contact discontinuity shows up, with more 
details appearing when increasing accuracy. 

This configuration leads to a complex flow structure, produced by a double 
Mach reflection of the shock at the wall. A self-similar flow develops as the 
fluid meets the reflecting wall. Two Mach stems develop, with two contact 
discontinuities. We have plotted in Fig. [13] the density contours at t = 0.2, 
when the main features have fully developed. The more challenging ones are 
the jet propagating to the right near the reflecting wall and the weak shock 
generated at the second Mach reflection, as seen in the enlarged area in Fig. [TH 

Our third-order results agree with the original ones [21] for the corresponding 
resolution: both the jet and the weak shock are clearly captured. Increasing 
both the resolution and the order-of-accuracy of the numerical algorithm, as 
shown in the subsequent panels in Fig. [TU we see more details of the jet 
roUing-up. Also, a vortex structure appears near the bottom wall, which starts 
affecting the diagonal contact discontinuity arising from the triple point. These 
high-resolution features, appearing in the last panel in Fig. [TH agree with the 
ones obtained with a WENO method of the same order (but double resolution, 
1/960) in Ref. [25]. This also agrees with the results of recent spectral (finite) 
volume simulations [26] , in which those structures show up gradually, as one is 
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getting more accurate simulations. This is another example in which a higher- 
order algorithm can be preferred, as it captures more detailed features of 
complex structures for a given resolution. 



7 Conclusions and outlook 

The numerical experiments presented in this paper provide clear evidence for 
a TVB behavior of the proposed schemes. This means that the total varia- 
tion growth is uniformly bounded, independently of the resolution, for a fixed 
evolution time. Moreover, the experimental pattern is a sudden growth of the 
total variation, which provides a time-independent bound for the rest of the 
evolution. This growth is confined to the mesh points directly connected with 
non-sonic critical points, especially near discontinuities. But the resulting rid- 
dles do not spread over smooth regions and the overall features of the solution 
are preserved as a result. In the case of compound shocks, however, the nu- 
merical simulations actually mystify the physical solution: the spurious riddles 
affect the contact point between the shock and the adjacent rarefaction wave, 
breaking the compound structure, even if the TVB behavior is still preserved. 

The proposed schemes are obtained from the unlimited version of the Osher- 
Chakrabarthy [6] linear flux-modification algorithms. The robustness of the 
unlimited version is related with the high compression factor of this algorithms 
family. We have actually improved the available estimates up to a remarkable 
value of 6 = 5, for the third-order case. This suggests that these estimates 
could be even improved by using alternative bound-setting procedures. Un- 
fortunately, even in the scalar case, we are not able to prove rigorously the 
TVB properties of these methods, although we are currently working in this 
direction. 

We have combined the unlimited Osher-Chakrabarthy algorithm with the sim- 
ple LLF flux formula. As a result, we have been able to derive the compact 
finite-difference scheme f l35p . which is equivalent to the corresponding finite- 
volume implementation in the unlimited case. This provides an extremely 
cost-efficient algorithm for dealing with the most common problems, even in 
presence of interacting dynamical shocks, as we have done in the Orszag-Tang 
2D vortex and the double Mach reflection cases. Of course, its use should be 
limited to convex-fiux problems, where compound shocks do not arise. 

The resulting finite-difference formula fl5^ is similar to the ones obtained 
by the 'artificial viscosity' approach (see for instance ref. [27]). The main 
difference is that the spectral radius plays a key role here in the dissipation 
term, providing some sort of 'adaptive viscosity'. Even in the most simple 
constant-speed case we can still see that the order of accuracy in our case 



23 



is dictated by the dissipation term, in contrast with the extra freedom one 
gets in the standard artificial viscosity approach. Moreover, our compression 
factor estimates provide specific prescriptions for the value of the dissipation 
coefficient. 
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